function res = predictiveVolRegression(Data,m,version)

% The egression reads:
% vol of excess bond returnt+m,k = alpha + (y_t_k-y_t_m) + eps_t+1
%
% INPUT: Must be a data matrix of dim T * numYields, 
% where yields have maturities of 1 quarter, 2 quarters,3 quarters, ...
[T,n]      = size(Data);
matSelect  = m:m:n;
shorRate   = Data(1:T-m,m);
yieldsUsed = Data(:,matSelect);
numYields  = size(yieldsUsed,2);
Betta      = NaN(2,numYields);
Betta_se   = NaN(2,numYields);
Betta_tstat= NaN(2,numYields);
Ydata      = NaN(T-m,numYields);
Xdata      = NaN(T-m,numYields);
R2         = NaN(1,numYields);
for i=1:numYields-1
   % Note that m*i gives k-m in relation to the notation above
   Y  = (-m*i/4*yieldsUsed(1+m:T,i) + m*(i+1)/4*yieldsUsed(1:T-m,i+1) - m/4*shorRate);
   x  = (yieldsUsed(1:T-m,i+1)-shorRate);
   X  = [ones(T-m,1) x];
   if version == 1
        res1 = runGARCHmodel(Y);
        volY = res1.conVol;
   elseif version == 2
        % GARCH or GARCH-AR model to measure uncertinaty in excess returns
        res1 = runGARCHmodel(Y);
        volY = res1.conVol;
        %res2 = runARGARCHmodel(Y,x,res1.params);
        %volY = res2.conVol;
   elseif version == 3
        % AR(1) model to measure uncertainty in excess returns
        tmp = nwest(Y(2:end),[ones(length(Y(2:end)),1) Y(1:end-1)],m);
        volY = [ 0; abs(tmp.resid)];
   elseif version == 4
       % GARCH or GARCH-AR model to measure uncertinaty in yields
        res1 = runGARCHmodel([yieldsUsed(1+m:T,i)-yieldsUsed(1+m-1:T-1,i)]);
        volY = res1.conVol;
        %res2 = runARGARCHmodel(yieldsUsed(1+m:T,i),NaN,res1.params);
        %volY = res2.conVol;
   elseif version == 5
       % AR(1) model to measure uncertainty in yields
       tmp  = nwest(yieldsUsed(1+m:T,i),[ones(T-m,1) yieldsUsed(1+m-1:T-1,i)],m);
       volY = abs(tmp.resid);
   end
   
   resOLS            = nwest(volY,X,m*2);
   Betta(:,i+1)      = resOLS.beta;
   Betta_tstat(:,i+1)= resOLS.tstat;
   Betta_se(:,i+1)   = resOLS.beta./resOLS.tstat;
   Ydata(:,i+1)      = volY;
   Xdata(:,i+1)      = x;
   R2(1,i+1)         = resOLS.rsqr;
end
% Output
res.betta       = Betta;
res.maturities  = matSelect;
res.betta_se    = Betta_se;
res.betta_tstat = Betta_tstat;
res.Ydata       = Ydata;
res.Xdata       = Xdata;
res.R2          = R2;

